Binary neutron stars: Equilibrium models beyond spatial conformal flatness 
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Equilibria of binary neutron stars in close circular orbits are computed numerically in a waveless 
formulation: The full Einstein-relativistic-Euler system is solved on an initial hypersurface to obtain 
an asymptotically flat form of the 4-metric and an extrinsic curvature whose time derivative vanishes 
in a comoving frame. Two independent numerical codes are developed, and solution sequences that 
model inspiraling binary neutron stars during the final several orbits are successfully computed. The 
binding energy of the system near its final orbit deviates from earlier results of third post-Newtonian 
and of spatially conformally flat calculations. The new solutions may serve as initial data for merger 
simulations and as members of quasiequilibrium sequences to generate gravitational wave templates, 
and may improve estimates of the gravitational- wave cutoff frequency set by the last inspiral orbit. 
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Introduction: Equilibria of close binary neutron stars 
in circular orbits, constructed numerically, have been 
studied as a model of the final several orbits of binary 
inspiral prior to merger (see [l| for a review). These nu- 
merical solutions have been used as initial data sets for 
merger simulations Q ; in quasi-equilibrium sequences, to 
estimate gravitational waveforms |3j,|4j; and to determine 
the cutoff frequency of the inspiral waves [f| . 

To maintain equilibrium circular orbits in general rel- 
ativity one must introduce an approximation that elim- 
inates the back reaction of gravitational radiation. An 
ansatz of this kind is the waveless approximation pro- 
posed by Isenberg |(| . One of his proposals was to choose 
a conformally flat spatial geometry maximally embedded 
in a spacetime. As a result the gravitational field is no 
longer dynamical; field equations for the metric compo- 
nents become elliptic equations. Wilson and Mathews 
later rediscovered this waveless approximation and ap- 
plied it to numerical computations of binary inspirals 7] . 

Although the Isenberg- Wilson-Mathews (IWM) for- 
mulation has been widely used for modeling binary neu- 
tron star and binary black hole inspiral in the past decade 
HI @- El EH 113 1 the error associated with its con- 
formally flat 3-geometry was studied only for stationary 
axisymmetric systems |l3| . In models of binary neutron 
stars, the estimated error in the orbital angular veloc- 
ity, Q,, is several percent 0, H3. implying a comparable 
deviation from circular orbits [15J. New waveless formu- 
lations, incorporating a generic form of the metric, are 
suitable for accurate computation of binary compact ob- 
jects 

Eiia. 

In this letter, we present the first results 
of numerical computations for binary neutron stars mod- 
eled in one of these formulations |l7j |. 

Formulation of the waveless spacetime: The new for- 
mulation 01 exactly solves the Einstein-Euler system 



written in 3+1 form on a s pac elike hypersurface. We 
follow notation |23j used in |l7j. The spacetime A4 = 
I x S is foliated by the family of spacelike hypersur- 
faces, E t = {£} x S. The future-pointing normal n a to 
Et is related to the timelike vector t a (the tangent dt to 
curves t — > (t, x), x G E) by t a = an a + (3 a , where a is 
the lapse, and where the shift j3 a satisfies (3 a n a = 0. A 
spatial metric j a b(t) defined on E t is equal to the projec- 
tion tensor j a p — g a p + n a np restricted to E t . In terms 
of a conformal factor ip and a conformally rescaled spa- 
tial metric j a b — 4>~ A lab, the metric g a p takes the form, 
ds 2 = -a 2 dt 2 + ip% 3 (dx l + ftdt)(dx j + ftdt), in a chart 
{t, x 1 }. A condition to specify the conformal decomposi- 
tion is det jab — det f a b, where f a b is a flat metric. 

In our waveless formulation we impose, as coordi- 
nate conditions, maximal slicing (K = 0) and the spa- 

o 

tially transverse condition Dbj = (the Dirac gauge 

[I7I Hi|). where Db is the covariant derivative with re- 
spect to the flat metric. We then restrict time-derivative 
terms in this gauge to guarantee that all components 
of the field equation are elliptic equations, and hence 
that all metric components, including the spatial met- 
ric, have Coulomb- type fall off While it is found 
to be sufficient to impose the condition dtj ab = 0(r~ 3 ) 
to have Coulomb-type fall off in the asymptotics, we im- 
pose a stronger condition: dtj ah — 0. For the other 
quantities, we impose helical symmetry: spacetime and 
fluid variables are dragged along by the helical vector 
k a = t a + il<f> a . For example, the time derivative of ex- 
trinsic curvature K a b is expressed as dfK a b = —£n<pK a b. 
The resulting field equations are solved on a slice Eo- 
The Hamiltonian constraint, momentum constraint, spa- 
tial trace and spatial tracefree part of the Einstein equa- 
tion are, respectively, regarded as elliptic equations for 
-0, /3 a , a and h a b := j a b — fab, while the extrinsic cur- 
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vature, K a },, for this formulation is computed from the 
metric components, K ab = j^Zpjab + ^7afc-D c (S10 c ). 

To compute the motion of binary neutron stars in cir- 
cular orbits, the flow field is assumed to be stationary 
in the rotating frame. Since any solution to the wave- 
less formulation satisfies all constraint equations, it is, 
in particular, an initial data set for the Einstein-Euler 
system. When one evolves such a binary neutron star 
solution by integrating the Einstein-Euler system, the 
orbits will deviate from exact circularity because of the 
radiation reaction force. Instead, one can construct an 
artificial spacetime with circular orbits by dragging the 
waveless solution on So along the vector k a = t a + fl(f> a , 
so that the spacetime has helical symmetry. Although 
the spacetime so constructed will not exactly satisfy Ein- 
stein's equation, a family of such spacetimes, associated 
with circular orbits of decreasing separation, will model 
the inspiral of a binary neutron star system during its fi- 
nal several orbits. Explicit forms of all equations for the 
fields and the matter are found in [ttI Il8j. 

Numerical methods: We have developed two indepen- 
dent numerical schemes to compute binary neutron star 
solutions. One is based on a finite difference method [Io| . 
the other one on a spectral method implemented via the 
C++ library Lorene Detailed convergence tests 

and calibration of each method will be published sepa- 
rately. In this letter, we show quantitative agreement 
of the two methods for h a b, which is the significant and 
reliable calibration for the new numerical solutions. 

In both methods, equations are written in Cartesian 
coordinate components, and they are solved numerically 
on spherical coordinate grids, r, 9, and 4>. In the finite 
difference method, an equally spaced grid is used from the 
center of orbital motion to 5i?o where there are n r = 16, 
24, and 32 grid points per Rq; from 5i?o to 10 4 i?o a 
logarithmically spaced grid has 60, 90, and 120 points 
(depending on the resolution). Here Rq is the geometric 
radius of a neutron star along a line passing through the 
center of orbit to the center of a star. Accordingly, for 9 
and 4> there are 32, 48, and 64 grid points each from 
to 7r/2 1(1 . For the spectral method, eight domains (a 
nucleus, six shells and a compactified domain extending 
up to infinity) around each star are used. In each domain, 
the number of collocation points is chosen to be N r xNgx 
= 25 x 17 x 16 and 33 x 21 x 20 0. 

Numerical solutions for binary neutron stars: A 
model of the evolutionary path of binary inspiral is given 
by a sequence of equilibria along which the neutron star 
matter is assumed to be isentropic; and the implied fluid 
flow is assumed to conserve the baryon number, entropy 
and vorticity of each fluid element Q, |2(| • In the case 
where the spins of component stars are negligible, the 
flow becomes irrotational; one can introduce the veloc- 
ity potential <I> by hu a — V Q $, where h is the specific 
enthalpy and u a is the fluid 4-velocity. For isentropic 
flow, one can assume a one-parameter equation of state, 




FIG. 1: Contours of (h xx — h yy )/2 in the a;y-plane, computed 
by the finite difference code (left) and by the spectral code 
(right). The binary separation a is given by a/Ro = 3.5. 
Contours extend from -0.014 to -0.002 with step 0.001. 
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FIG. 2: Components hij along the a;-axis, normalized by 
X — tv/Q. A neutron star extends from x/X = 0.02024 to 
0.07422. Curves labeled FD and SP display results of the 
finite difference and spectral codes, respectively. 



V = v{p)i w ith p the baryon mass density. The matter is 
then described by two independent variables, a thermo- 
dynamic variable such as p/p, and the velocity potential 
$. In this letter, we assume a polytropic equation of 
state p = Kp T with adiabatic index r = 2, and we dis- 
play results for equal-mass binaries with the rest mass 
of each star to be that of a single spherical star of com- 
pactness (M/R)^ = 0.17. (Note: The maximum com- 
pactness of a spherical star for this equation of state is 
(M/Rjoa = 0.216. The compactness (M/R)^ is defined 
as the ratio of graviational mass to circumferential radius 
of an isolated spherical star with the same rest mass.) 

In Fig. ^ contours of the components hij computed by 
the two numerical codes are shown for selected solutions. 
In these solutions, the separation in coordinate distance 
between the coordinate center of each neutron star is set 
to a/Ro — 3.5. From these contours, one can verify qual- 
itative agreement of the results from the two independent 
numerical methods. In Fig. EJ components hij along the 
x-axis are plotted for the same solution, where the cc-axis 
passes through the centers of the neutron stars. Preci- 
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FIG. 3: Virial error vs. angular velocity Q, normalized by 
Moo , twice the gravitational mass of an isolated neutron star. 
Each curve labeled FD shows results of a finite difference code 
with a given resolution. Curves labeled SP and IWM show 
results of the spectral code and the spatially conformally flat 
approximation, respectively. 



sion of integral quantities characterizing the solutions is 
shown by the finite difference (spectral) method compar- 
isons: flMoo = 0.03565 (0.03565), M AD m/Moo = 0.98825 
(0.98826). and J/M^ = 0.9212 (0.9165). 

In ^3] j h is shown that the ADM mass, Madm, and the 
asymptotic Komar mass, Mk are equal, Madm = Mk, 
under asymptotic conditions satisfied by the solutions in 
the present formulation. The equality is related to a virial 
relation for the equilibrium, 



I); 



fVpTj^f^cPx = 0, 



(1) 



that we use to evaluate the accuracy of numerical so- 
lutions. Fig. |21 shows the computed value of the virial 
integral in Eq. lfT)l. normalized by Madm, along the se- 
quence. We also evaluated Madm and Mk each defined 
by the surface integral in the asymptotics, and confirmed 
that, for each model, the difference of the two masses is 
no larger than |Madm — Mk|/Madm ~ 0.01% for the 
finite difference method and ~ 0.001% for the spectral 
method; these errors are consistent with the numerical 
errors of the virial relation shown in Fig. |21 

In Fig. 01 the binding energy E b = Madm — M M along 
the sequence is plotted and compared with that resulting 
from a third post-Newtonian (3PN) calculation [2l| and 
IWM formulation. The waveless sequence fits the 3PN 
curve well at larger separation, and reaches a configura- 
tion with a cusp without any turning point in the binding 
energy curve, in agreement with results of the IWM for- 
mulation [3, Uij (the spectral code does not yet converge 
for the closest orbits - largest QM^- of Figs. |21 and 0] 
because it is more sensitive to tidal deformation: higher 
multipoles in the density of each star lead to a divergent 
iteration). The binding energy E b of the waveless se- 
quences clearly deviates from that of the 3PN and IWM 
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FIG. 4: Binding energy Eb := Madm — Moo with respect to 
the normalized angular velocity along the sequence. Curves 
are labeled as in Fig. 01 The thin solid curve corresponds to 
the third order post-Newtonian calculation l2l| . 
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FIG. 5: Norm of Lie derivative of the extrinsic curvature along 
the helical vector [cf. Eq. JJ] with respect to the normalized 
angular velocity. 



sequences at the larger values of OMqo. This suggests 
that the 3PN and IWM formulations each overestimate 
the binding energy - in the 3PN case, by neglecting the 
tidal deformation, in the IWM formulation by neglecting 
the contribution from h ab . 

Finally, to estimate the deviation of the orbit from cir- 
cularity, we evaluate the formal expression for the extrin- 
sic curvature of a solution with exact helical symmetry, 
(the case for which the time-evolved data has an exactly 
circular orbit), K ab = £ p +0,4,1 ab- Because £kK ab van- 
ishes for exact helical symmetry, its norm, defined on the 
support V of the fluid, 



£k.Kn 



r c r d £kK ab £ k K cd ^d 6 



1/2 



(2) 



is a measure of the deviation from circularity. 

Fig. El shows that, for all separations, the values of 
II £kK ab || for the waveless solutions are more than an 
order of magnitude smaller than those of IWM solutions. 
The result supports the expectation that IWM data en- 
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forces circularity with significantly less accuracy than the 
corresponding waveless solutions, even for larger separa- 
tion. This may be important: Even for a sudden turn- 
on of radiation-reaction a post-Newtonian analysis [l5| 
shows eccentricity < 1.5% at Q,M < 0.03 for initially 
circular orbits, and the gravitational radiation reaction 
should be more gradual for our waveless data sets. 

Discussion: In 2nd post-Newtonian theory (e.g. |l4j|L 
the correction AEb to the binding energy due to the 
contribution of h a b is of order M 00 h a bV a v b , where the 
magnitude of orbital velocity, v a , may be typically v « 
0.34(fiM oo /0.04) 1 / 3 . Since h ab is 0(v 4 ), AEb/M^ = 
0{v 6 ) ~ 10~ 3 for CWoo - 0.04. This agrees with the 
difference between the binding energies calculated by the 
IWM and waveless formulations in Fig. 0] 

The quantity dEb /dfl is important for the data analysis 
of gravitational waves, because it determines the evolu- 
tion of the wave's phase, $gw = 2 J fl(t)dt. In adiabatic 
evolution, the time dependence of angular velocity f2(t) is 
calculated from dfl/dt = \(dE / dt) gw\/ {dEb /dd), where 
(dE/dt)Gw is the luminosity of gravitational waves. Our 
present result shows that the derivative dEb /dQ, of wave- 
less sequences is ~ 10-15% larger than those of IWM 
and 3PN curves for QMoa > 0.035. Since ~ 2 orbits 
are maintained from QMoq = 0.035 to merger for the 
case with (M/R)^ = 0.17 Q, the error in the IWM and 
3PN values of $gw would accumulate to ~ 50% over the 
last ~ 2 orbits. The phase error leads to error during 
the final orbits before merger of the computed frequency, 
whose final behavior constrains the equation of state of 
nuclear matter j^H^]. Waveless solutions may determine 
phase and frequency with significantly greater accuracy - 
particularly if, to overcome radial-motion error, one first 
calibrates the frequencies of a set of quasiequililbrium 
sequences, using (for example) time evolutions. 

Phase error may be much larger for the final orbits of 
binary black hole and black hole-neutron star inspirals. 
In these cases, QM^ in the last orbit may reach or exceed 
0.1 (e.g. HHEl). Since AE b is of order 0(t> 6 ), the phase 
error is likely to be of order unity for QMoo > 0.1. There- 
fore, a template constructed from the IWM formulation 
may cause a systematic error in the data analysis. Our 
waveless approximation may improve binary black hole 
and black hole-neutron star solutions for this purpose. 
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